home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / POLYROOT.ZIP / PRBM.FOR < prev    next >
Text File  |  1985-11-29  |  8KB  |  248 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE PRBM
  5. C
  6. C        PURPOSE
  7. C           TO CALCULATE ALL REAL AND COMPLEX ROOTS OF A GIVEN
  8. C           POLYNOMIAL WITH REAL COEFFICIENTS.
  9. C
  10. C        USAGE
  11. C           CALL PRBM (C,IC,RR,RC,POL,IR,IER)
  12. C
  13. C        DESCRIPTION OF PARAMETERS
  14. C           C      - INPUT VECTOR CONTAINING THE COEFFICIENTS OF THE
  15. C                    GIVEN POLYNOMIAL. COEFFICIENTS ARE ORDERED FROM
  16. C                    LOW TO HIGH. ON RETURN COEFFICIENTS ARE DIVIDED
  17. C                    BY THE LAST NONZERO TERM.
  18. C           IC     - DIMENSION OF VECTORS C, RR, RC, AND POL.
  19. C           RR     - RESULTANT VECTOR OF REAL PARTS OF THE ROOTS.
  20. C           RC     - RESULTANT VECTOR OF COMPLEX PARTS OF THE ROOTS.
  21. C           POL    - RESULTANT VECTOR OF COEFFICIENTS OF THE POLYNOMIAL
  22. C                    WITH CALCULATED ROOTS. COEFFICIENTS ARE ORDERED
  23. C                    FROM LOW TO HIGH (SEE REMARK 4).
  24. C           IR     - OUTPUT VALUE SPECIFYING THE NUMBER OF CALCULATED
  25. C                    ROOTS. NORMALLY IR IS EQUAL TO IC-1.
  26. C           IER    - RESULTANT ERROR PARAMETER CODED AS FOLLOWS
  27. C                     IER=0  - NO ERROR,
  28. C                     IER=1  - SUBROUTINE PQFB RECORDS POOR CONVERGENCE
  29. C                              AT SOME QUADRATIC FACTORIZATION WITHIN
  30. C                              50 ITERATION STEPS,
  31. C                     IER=2  - POLYNOMIAL IS DEGENERATE, I.E. ZERO OR
  32. C                              CONSTANT,
  33. C                              OR OVERFLOW IN NORMALIZATION OF GIVEN
  34. C                              POLYNOMIAL,
  35. C                     IER=3  - THE SUBROUTINE IS BYPASSED DUE TO
  36. C                              SUCCESSIVE ZERO DIVISORS OR OVERFLOWS
  37. C                              IN QUADRATIC FACTORIZATION OR DUE TO
  38. C                              COMPLETELY UNSATISFACTORY ACCURACY,
  39. C                     IER=-1 - CALCULATED COEFFICIENT VECTOR HAS LESS
  40. C                              THAN THREE CORRECT SIGNIFICANT DIGITS.
  41. C                              THIS REVEALS POOR ACCURACY OF CALCULATED
  42. C                              ROOTS.
  43. C
  44. C        REMARKS
  45. C           (1) REAL PARTS OF THE ROOTS ARE STORED IN RR(1) UP TO RR(IR)
  46. C               AND CORRESPONDING COMPLEX PARTS IN RC(1) UP TO RC(IR).
  47. C           (2) ERROR MESSAGE IER=1 INDICATES POOR CONVERGENCE WITHIN
  48. C               50 ITERATION STEPS AT SOME QUADRQTIC FACTORIZATION
  49. C               PERFORMED BY SUBROUTINE PQFB.
  50. C           (3) NO ACTION BESIDES ERROR MESSAGE IER=2 IN CASE OF A ZERO
  51. C               OR CONSTANT POLYNOMIAL. THE SAME ERROR MESSAGE IS GIVEN
  52. C               IN CASE OF AN OVERFLOW IN NORMALIZATION OF GIVEN
  53. C               POLYNOMIAL.
  54. C           (4) ERROR MESSAGE IER=3 INDICATES SUCCESSIVE ZERO DIVISORS
  55. C               OR OVERFLOWS OR COMPLETELY UNSATISFACTORY ACCURACY AT
  56. C               ANY QUADRATIC FACTORIZATION PERFORMED BY
  57. C               SUBROUTINE PQFB. IN THIS CASE CALCULATION IS BYPASSED.
  58. C               IR RECORDS THE NUMBER OF CALCULATED ROOTS.
  59. C               POL(1),...,POL(J-IR) ARE THE COEFFICIENTS OF THE
  60. C               REMAINING POLYNOMIAL, WHERE J IS THE ACTUAL NUMBER OF
  61. C               COEFFICIENTS IN VECTOR C (NORMALLY J=IC).
  62. C           (5) IF CALCULATED COEFFICIENT VECTOR HAS LESS THAN THREE
  63. C               CORRECT SIGNIFICANT DIGITS THOUGH ALL QUADRATIC
  64. C               FACTORIZATIONS SHOWED SATISFACTORY ACCURACY, THE ERROR
  65. C               MESSAGE IER=-1 IS GIVEN.
  66. C           (6) THE FINAL COMPARISON BETWEEN GIVEN AND CALCULATED
  67. C               COEFFICIENT VECTOR IS PERFORMED ONLY IF ALL ROOTS HAVE
  68. C               BEEN CALCULATED. IN THIS CASE THE NUMBER OF ROOTS IR IS
  69. C               EQUAL TO THE ACTUAL DEGREE OF THE POLYNOMIAL (NORMALLY
  70. C               IR=IC-1). THE MAXIMAL RELATIVE ERROR OF THE COEFFICIENT
  71. C               VECTOR IS RECORDED IN RR(IR+1).
  72. C
  73. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  74. C           SUBROUTINE PQFB     QUADRATIC FACTORIZATION OF A POLYNOMIAL
  75. C                               BY BAIRSTOW ITERATION.
  76. C
  77. C        METHOD
  78. C           THE ROOTS OF THE POLYNOMIAL ARE CALCULATED BY MEANS OF
  79. C           SUCCESSIVE QUADRATIC FACTORIZATION PERFORMED BY BAIRSTOW
  80. C           ITERATION. X**2 IS USED AS INITIAL GUESS FOR THE FIRST
  81. C           QUADRATIC FACTOR, AND FURTHER EACH CALCULATED QUADRATIC
  82. C           FACTOR IS USED AS INITIAL GUESS FOR THE NEXT ONE. AFTER
  83. C           COMPUTATION OF ALL ROOTS THE COEFFICIENT VECTOR IS
  84. C           CALCULATED AND COMPARED WITH THE GIVEN ONE.
  85. C           FOR REFERENCE, SEE J. H. WILKINSON, THE EVALUATION OF THE
  86. C           ZEROS OF ILL-CONDITIONED POLYNOMIALS (PART ONE AND TWO),
  87. C           NUMERISCHE MATHEMATIK, VOL.1 (1959), PP.150-180.
  88. C
  89. C     ..................................................................
  90. C
  91.       SUBROUTINE PRBM(C,IC,RR,RC,POL,IR,IER)
  92. C
  93. C
  94.       DIMENSION C(1),RR(1),RC(1),POL(1),Q(4)
  95. C
  96. C        TEST ON LEADING ZERO COEFFICIENTS
  97.       EPS=1.E-3
  98.       LIM=50
  99.       IR=IC+1
  100.     1 IR=IR-1
  101.       IF(IR-1)42,42,2
  102.     2 IF(C(IR))3,1,3
  103. C
  104. C        WORK UP ZERO ROOTS AND NORMALIZE REMAINING POLYNOMIAL
  105.     3 IER=0
  106.       J=IR
  107.       L=0
  108.       A=C(IR)
  109.       DO 8 I=1,IR
  110.       IF(L)4,4,7
  111.     4 IF(C(I))6,5,6
  112.     5 RR(I)=0.
  113.       RC(I)=0.
  114.       POL(J)=0.
  115.       J=J-1
  116.       GO TO 8
  117.     6 L=1
  118.       IST=I
  119.       J=0
  120.     7 J=J+1
  121.       C(I)=C(I)/A
  122.       POL(J)=C(I)
  123.       CALL OVERFL(N)
  124.       IF(N-2)42,8,8
  125.     8 CONTINUE
  126. C
  127. C        START BAIRSTOW ITERATION
  128.       Q1=0.
  129.       Q2=0.
  130.     9 IF(J-2)33,10,14
  131. C
  132. C        DEGREE OF RESTPOLYNOMIAL IS EQUAL TO ONE
  133.    10 A=POL(1)
  134.       RR(IST)=-A
  135.       RC(IST)=0.
  136.       IR=IR-1
  137.       Q2=0.
  138.       IF(IR-1)13,13,11
  139.    11 DO 12 I=2,IR
  140.       Q1=Q2
  141.       Q2=POL(I+1)
  142.    12 POL(I)=A*Q2+Q1
  143.    13 POL(IR+1)=A+Q2
  144.       GO TO 34
  145. C        THIS IS BRANCH TO COMPARISON OF COEFFICIENT VECTORS C AND POL
  146. C
  147. C        DEGREE OF RESTPOLYNOMIAL IS GREATER THAN ONE
  148.    14 DO 22 L=1,10
  149.       N=1
  150.    15 Q(1)=Q1
  151.       Q(2)=Q2
  152.       CALL PQFB(POL,J,Q,LIM,I)
  153.       IF(I)16,24,23
  154.    16 IF(Q1)18,17,18
  155.    17 IF(Q2)18,21,18
  156.    18 GO TO (19,20,19,21),N
  157.    19 Q1=-Q1
  158.       N=N+1
  159.       GO TO 15
  160.    20 Q2=-Q2
  161.       N=N+1
  162.       GO TO 15
  163.    21 Q1=1.+Q1
  164.    22 Q2=1.-Q2
  165. C
  166. C        ERROR EXIT DUE TO UNSATISFACTORY RESULTS OF FACTORIZATION
  167.       IER=3
  168.       IR=IR-J
  169.       RETURN
  170. C
  171. C        WORK UP RESULTS OF QUADRATIC FACTORIZATION
  172.    23 IER=1
  173.    24 Q1=Q(1)
  174.       Q2=Q(2)
  175. C
  176. C        PERFORM DIVISION OF FACTORIZED POLYNOMIAL BY QUADRATIC FACTOR
  177.       B=0.
  178.       A=0.
  179.       I=J
  180.    25 H=-Q1*B-Q2*A+POL(I)
  181.       POL(I)=B
  182.       B=A
  183.       A=H
  184.       I=I-1
  185.       IF(I-2)26,26,25
  186.    26 POL(2)=B
  187.       POL(1)=A
  188. C
  189. C        MULTIPLY POLYNOMIAL WITH CALCULATED ROOTS BY QUADRATIC FACTOR
  190.       L=IR-1
  191.       IF(J-L)27,27,29
  192.    27 DO 28 I=J,L
  193.    28 POL(I-1)=POL(I-1)+POL(I)*Q2+POL(I+1)*Q1
  194.    29 POL(L)=POL(L)+POL(L+1)*Q2+Q1
  195.       POL(IR)=POL(IR)+Q2
  196. C
  197. C        CALCULATE ROOT-PAIR FROM QUADRATIC FACTOR X*X+Q2*X+Q1
  198.       H=-.5*Q2
  199.       A=H*H-Q1
  200.       B=SQRT(ABS(A))
  201.       IF(A)30,30,31
  202.    30 RR(IST)=H
  203.       RC(IST)=B
  204.       IST=IST+1
  205.       RR(IST)=H
  206.       RC(IST)=-B
  207.       GO TO 32
  208.    31 B=H+SIGN(B,H)
  209.       RR(IST)=Q1/B
  210.       RC(IST)=0.
  211.       IST=IST+1
  212.       RR(IST)=B
  213.       RC(IST)=0.
  214.    32 IST=IST+1
  215.       J=J-2
  216.       GO TO 9
  217. C
  218. C        SHIFT BACK ELEMENTS OF POL BY 1 AND COMPARE VECTORS POL AND C
  219.    33 IR=IR-1
  220.    34 A=0.
  221.       DO 38 I=1,IR
  222.       Q1=C(I)
  223.       Q2=POL(I+1)
  224.       POL(I)=Q2
  225.       IF(Q1)35,36,35
  226.    35 Q2=(Q1-Q2)/Q1
  227.    36 Q2=ABS(Q2)
  228.       IF(Q2-A)38,38,37
  229.    37 A=Q2
  230.    38 CONTINUE
  231.       I=IR+1
  232.       POL(I)=1.
  233.       RR(I)=A
  234.       RC(I)=0.
  235.       IF(IER)39,39,41
  236.    39 IF(A-EPS)41,41,40
  237. C
  238. C        WARNING DUE TO POOR ACCURACY OF CALCULATED COEFFICIENT VECTOR
  239.    40 IER=-1
  240.    41 RETURN
  241. C
  242. C        ERROR EXIT DUE TO DEGENERATE POLYNOMIAL OR OVERFLOW IN
  243. C        NORMALIZATION
  244.    42 IER=2
  245.       IR=0
  246.       RETURN
  247.       END
  248.